Análisis 3
QS3
Nos vamos a quedar con la onda 3 y vamos a analizar la zona entre 65º y 40ºS y entre 925 y 30 hPa
g_data <- qs_gh_interp[QS == 3]
region_3 <- c(latmin = -65, latmax = -40, levmin = 30, levmax = 925)
M <- ceiling(g_data[, max(Amplitud, na.rm = T)])
lineas_width <- 10
lineas <- seq(0, M, by = lineas_width)
s <- 0.7
ggplot(g_data[modelo == "sp"], aes(lat, lev)) +
scale_y_continuous(trans = "reverselog") +
scale_x_reverse() +
facet_grid(month ~ .) +
geom_raster(data = g_data[modelo == "nc"],
aes(fill = cut(Amplitud, breaks = lineas))) +
geom_contour(aes(z = Amplitud), binwidth = lineas_width,
color = "black", size = s, alpha = 0.7) +
# geom_dl(aes(z = Amplitud, label = ..level..),
# binwidth = lineas_width, stat = "contour", method = "top.pieces") +
annotate(geom = "rect",
xmin = region_3[1], xmax = region_3[2],
ymin = region_3[3], ymax = region_3[4],
fill = NA, color = "black") +
scale_fill_brewer(palette = "Blues", name = "NCEP") +
theme_elio + ggtitle(paste0("Amplitud de QS3 por latitud y mes \nNCEP (sombreado) y SPEEDY (contornos) - lineas = ", lineas_width))
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'key.row'
## Note: no visible binding for global variable 'key.col'
## Note: no visible binding for global variable 'label.row'
## Note: no visible binding for global variable 'label.col'
## Note: no visible binding for global variable 'key.row'
## Note: no visible binding for global variable 'key.col'
## Note: no visible binding for global variable 'label.row'
## Note: no visible binding for global variable 'label.col'
## Note: no visible binding for global variable 'key.row'
## Note: no visible binding for global variable 'key.col'
## Note: no visible binding for global variable 'label.row'
## Note: no visible binding for global variable 'label.col'
## Note: no visible binding for global variable 'key.row'
## Note: no visible binding for global variable 'key.col'
## Note: no visible binding for global variable 'label.row'
## Note: no visible binding for global variable 'label.col'En esa caja puedo tomar el valor promedio (que sería equivalente a la integral, ya que el área se mantiene igual) o el máximo. El ciclo anual de ambas variables tanto para NCEP como para SPEEDY se muestra en la siguiente figura:
qs3_gh <- qs_gh[QS == 3, ]
invisible(
qs3_gh[, region := (lat >= region_3[1] & lat <= region_3[2] &
lev >= region_3[3] & lev <= region_3[4])]
)
g_data <- qs3_gh[region == T,
.(Mean_Ampl = mean(Amplitud), Max_Ampl = max(Amplitud)),
by = .(month, modelo)]
invisible(
g_data[, Area_equiv := Mean_Ampl/Max_Ampl]
)
g <- ggplot(g_data, aes(as.numeric(month), Area_equiv)) +
geom_line(aes(color = modelo)) +
scale_color_brewer(type = "qual", palette = 6, labels = c("NCEP", "SPEEDY"),
name = "Modelo", direction = -1) +
scale_x_continuous(breaks = 1:12, labels = month.abb, name = "Mes") +
scale_y_continuous(limits = c(NA, 1)) +
ylab("Amplitud media / Amplitud máxima") +
theme_elio
g_data <- melt(g_data[, !"Area_equiv", with = F],
id.vars = c("modelo", "month"),
variable.name = "Tipo", value.name = "Amplitud")
ggplot(g_data, aes(as.numeric(month), Amplitud)) +
geom_line(aes(color = modelo, linetype = Tipo)) +
scale_x_continuous(breaks = 1:12, labels = month.abb, name = "Mes") +
scale_color_brewer(type = "qual", palette = 6, labels = mod_labels,
name = "Modelo", direction = -1) +
theme_elio
Se observa que no hay mucha correlación entre los modelos. Mientras que NCEP tiene un ciclo semianual con máximos en febrero y agosto, Speedy tiene los máximos en junio y noviembre. Casi una anticorrelación perfecta.
Otra medida de posible relevancia es la relación entre la amplitud máxima y la media. Ésto da una idea de cuan concentrado está la anomalía de geopotencial. Valores cercanos a 1 implican que la misma está distribuida de forma pareja en todo el área de estudio, mientras que valores menores implican una distribución más concentrada.
g
Esta variable muestra un comportamiento muy parecido en ambos modelos, con un máximo en el invierno y un mínimo en verano, indicando que el patrón de QS3 se encuentra mas localizado en verano que en invierno.
Una cosa que hay que tener cuidado es cómo caracterizar la amplitud de la onda 3 en toda la región mediante un número para hacer una serie temporal. La elección más naive es la media, pero ésto sólo es útil si la variable tiene una distribución simétrica y unimodal. Una variable como la amplitud está acotado inferiormente por el cero, por lo que la primera propiedad no se cumple a priori (aunque puede ser aproximada si la amplitud es grande).
La funciones de distribución para cada mes estimadas a partir de los datos muestran que en varios meses tampoco se cumple la segunda propiedad.
ggplot(qs3_gh[region == T, ], aes(Amplitud)) +
geom_density(aes(color = modelo)) + facet_grid(month~.) +
scale_color_brewer(type = "qual", palette = 6, labels = mod_labels,
name = "Modelo", direction = -1) +
xlab("Densidad") +
theme_elio
## Note: no visible binding for global variable 'y'
## Note: no visible binding for global variable 'x'
## Note: no visible binding for global variable 'x'
## Note: no visible binding for global variable 'ymax'
## Note: no visible binding for global variable 'ymin'En particular, agosto y noviembre muestran distribuciones multimodales en speedy y ncep respectivamente. Esto quizás podría poner en duda la validez de usar la media como representativo de la actividad de la onda 3 en toda la región.
Serie temporal onda 3
Independientemente de esto, armo una serie temporal mensual con la amplitud de la QS3 para el período 1985-2015 usando datos de NCEP.
# load("NCEP/todo_ncep.rda")
#
# ncep_r3 <- ncep[lat >= region_3[1] & lat <= region_3[2] &
# lev >= region_3[3] & lev <= region_3[4]]
load("NCEP/ncep_r3.rda")
invisible({
ncep_r3[, agh := anom(gh), by = .(lat, lev, month)]
qs3 <- ncep_r3[, qs_fit(agh, n = 3), by = .(lat, lev, time)]
qs3[, month := as.numeric(stringi::stri_sub(time, 6, 7))]
qs3[, month := factor(month, levels = c(12, 1:11), ordered = T)]
qs3[, year := as.numeric(stringi::stri_sub(time, 1, 4))]
})Meses de verano donde activo e inactivo se definen si la amplitud es mayor, menor o que la media (del mes) en 1 desvío estandar (del mes):
qs3_mean <- qs3[, .(Mean_Ampl = mean(Amplitud),
Max_Ampl = max(Amplitud)),
by = .(month, year, time)]
invisible(
qs3_mean[, `:=`(MeanMean_Ampl = mean(Mean_Ampl),
SDMean_Ampl = sd(Mean_Ampl),
MeanMax_Ampl = mean(Max_Ampl),
SDMax_Ampl = sd(Max_Ampl)),
by = .(month)]
)
shift_season <- function(month, year) {
ifelse(month == 12, year + 1, year)
}
invisible({
qs3_mean[, Año_estacion := shift_season(month, year)]
qs3_mean[, Anom := ifelse(Mean_Ampl > MeanMean_Ampl + SDMean_Ampl,
"Activo",
ifelse(Mean_Ampl < MeanMean_Ampl - SDMean_Ampl,
"Inactivo", "Normal"))]
})
# Alternativa: usando terciles.
# qs3_mean[, Anom := cut(Mean_Ampl,
# breaks = quantile(Mean_Ampl, seq(0, 1, length.out = 4)),
# include.lowest = T, labels = c("Inactivo", "Normal", "Activo")),
# by = month]
years <- seq(1984, 2016, by = 1)
lab <- stri_sub(years, 3, 4)
labs <- paste0(lab, "/", shift(lab, type = "lead"))
n <- length(years)
years <- years[-c(1)]
labs <- labs[-c(n)]
# years <- unique(qs3_mean$Año_estacion)
g <- ggplot(qs3_mean,
aes(Año_estacion, Mean_Ampl)) +
geom_col(aes(fill = Anom)) +
geom_line(aes(y = MeanMean_Ampl), linetype = 2, color = "gray4") +
scale_x_continuous(breaks = jumpby(years, 2),
labels = jumpby(labs, 2),
name = "Año de fin de estación") +
scale_fill_brewer(palette = "Set1") +
ylab("Amplitud media QS3") +
theme_elio +
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1))
g + facet_grid(month~., labeller = labeller(month = month.abb))
## Note: no visible binding for global variable 'y'
## Note: no visible binding for global variable 'y'
## Note: no visible binding for global variable 'x'
## Note: no visible binding for global variable 'width'
## Note: no visible binding for global variable 'x'
## Note: no visible binding for global variable 'width'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'key.row'
## Note: no visible binding for global variable 'key.col'
## Note: no visible binding for global variable 'label.row'
## Note: no visible binding for global variable 'label.col'
## Note: no visible binding for global variable 'key.row'
## Note: no visible binding for global variable 'key.col'
## Note: no visible binding for global variable 'label.row'
## Note: no visible binding for global variable 'label.col'
## Note: no visible binding for global variable 'key.row'
## Note: no visible binding for global variable 'key.col'
## Note: no visible binding for global variable 'label.row'
## Note: no visible binding for global variable 'label.col'
## Note: no visible binding for global variable 'key.row'
## Note: no visible binding for global variable 'key.col'
## Note: no visible binding for global variable 'label.row'
## Note: no visible binding for global variable 'label.col'La gran mayoría de los meses de veranos son normales con 5-6 veranos activo o inactivos.
Se observa que no hay muchos meses con amplitud mucho mayor que la media y que no están agrupados en un mismo verano. Es decir, no hay mucha correlación entre la amplitud de un mes y el siguiente.
Para ver esto último con más exactitud, se puede estimar la función de autocorrelación.
autocor <- acf.sig(qs3_mean$Mean_Ampl, method = "large.lag")
ggplot(autocor, aes(lag, acf)) +
geom_line() + geom_point() +
geom_line(aes(y = sig.cut), linetype = 2) +
geom_line(aes(y = -sig.cut), linetype = 2) +
geom_hline(yintercept = 0) +
scale_x_continuous(limits = c(0, 12*3),
name = "Lag (meses)",
breaks = seq(0, 12*3, by = 6)) +
ylab("Autocorrelación mensual QS3") +
theme_elio
## Note: no visible binding for global variable 'xend'
## Note: no visible binding for global variable 'yend'
## Note: no visible binding for global variable 'x'
## Note: no visible binding for global variable 'y'Se ve que la autocorrelación a lag 1 no es significativa (utilizando método large lag), pero además se ve un evidente ciclo anual.
ggplot(qs3_mean, aes(factor(month, labels = month.abb), Mean_Ampl)) +
geom_boxplot() +
xlab("Mes") + ylab("Amplitud media QS3") +
theme_elio
## Note: no visible binding for global variable 'weight'
## Note: no visible binding for global variable 'xmin'
## Note: no visible binding for global variable 'xmax'
## Note: no visible binding for global variable 'y'
## Note: no visible binding for global variable 'size'¿Tiene otros ciclos? Le quito el ciclo anual restando la media anual y hago fourier.
fourier <- as.dt(convert.fft(fft(qs3_mean[, .(Mean_Ampl - MeanMean_Ampl)]$V1)))
ggplot(fourier, aes(freq, ampl)) +
geom_line() +
theme_elio + xlab("Frecuencia") + ylab("Amplitud")
No, no hay ningún ciclo obvio.
qs3 <- merge(qs3, qs3_mean[, .(month, year, Anom)])
#
# g_data <- qs3[, .(Amplitud = mean(Amplitud)), by = .(lat, lev, Anom, month)]
#
# g_data <- g_data[, interp.dt(lat, lev, Amplitud, linear = T,
# yo = interp.levs),
# by = .(Anom, month)]
# colnames(g_data) <- c("Anom", "month", "lat", "lev", "Amplitud")
# # g_data <- g_data[month %in% 1:2]
# ggplot(g_data, aes(lat, lev)) +
# scale_y_continuous(trans = "reverselog") +
# geom_contour(aes(z = Amplitud), binwidth = 10,
# color = "black") +
# geom_dl(aes(z = Amplitud, label = ..level..),
# stat = "contour", binwidth = 20, color = "black",
# method = "top.pieces") +
# facet_grid(month~Anom, labeller = labeller(month = month.abb)) +
# theme_elio
# Esto es la distribución espacial de la amplitud en la región de estudio para cada mes y clasificación. No sé si sirve para mucho. Se ve que en febrero en algunos meses hay más separación que en otros.
ggplot(qs3, aes(Amplitud)) +
geom_density(aes(color = Anom)) +
facet_wrap(~month, ncol = 4, labeller = labeller(month = month.abb)) +
theme_elio
## Note: no visible binding for global variable 'y'
## Note: no visible binding for global variable 'x'
## Note: no visible binding for global variable 'x'
## Note: no visible binding for global variable 'ymax'
## Note: no visible binding for global variable 'ymin'Composiciones de variables.
Ahora la idea es ver los campos de variables que podrían ser forantes y comparar años activos vs. inactivos.
SST
load("NCEP/sst_todo.rda")
## Agrego máscara de océanos
library(maptools)
library(maps)
make_mask <- function(lat, lon) {
seamask <- map("world2", fill=TRUE, col = "transparent", plot = F)
IDs <- sapply(strsplit(seamask$names, ":"), function(x) x[1])
seamask <- map2SpatialPolygons(seamask, IDs = IDs,
proj4string = CRS("+proj=longlat +datum=WGS84"))
points <- SpatialPoints(expand.grid(lon, lat),
proj4string = CRS(proj4string(seamask)))
sea <- is.na(over(points, seamask))
return(sea)
}
lat <- unique(sst$lat)
lon <- unique(sst$lon)
sea <- make_mask(lat, lon)
invisible(sst[, sea := sea])
sst <- merge(sst, qs3_mean[, .(time, Anom, Mean_Ampl)], all = T)
invisible({
sst[, month := as.numeric(stringi::stri_sub(time, 6, 7))]
sst[, month := factor(month, levels = c(12, 1:11), ordered = T)]
sst[sea == F, sst := NA]
sst[, msst := mean(sst, na.rm = T), by = .(lat, month)]
sst[, asst := sst - msst]
})Probé la composición de SST de años activos e inactivos y prácticamente no se ve diferencia, así que directamente muestro la diferencia (activo - inactivo):
sst_comp <- sst[!is.na(Anom),
.(sst = mean(sst, na.rm = T)),
by = .(lat, lon, Anom, month)]
tmp <- sst_comp[, .(Anom = "Diferencia",
sst = sst[Anom == "Activo"] - sst[Anom == "Inactivo"]),
by = .(lat, lon, month)]
sst_comp <- rbind(sst_comp, tmp)
remove(tmp)
gdata <- sst_comp[Anom == "Diferencia"]
ggplot(gdata, aes(lon, lat)) +
geom_contour(aes(z = sst, color = ..level..), binwidth = .5) +
world3 +
coord_quickmap() +
scale_color_gradient2(high = muted("red"), low = muted("blue"), name = "SST \n(Activo - Inactivo)") +
theme_elio + facet_wrap(~month, labeller = labeller(month = month.abb),
ncol = 3) +
labs(title = "Diferencia de SST entre años activos e inactivos",
subtitle = "lineas = 0.5")
invisible(
sst[!is.na(Mean_Ampl) & !is.na(asst), sd := sd(asst), by = .(lat, lon, month)]
)
# sst_lm <-
# sst[!is.na(Mean_Ampl) & !is.na(asst) & sd != 0,
# {
# a <- summary(lm(Mean_Ampl ~ asst))
# list(estimate = a$coefficients[2, 1],
# se = a$coefficients[2, 2])
# },
# by = .(lon, lat, month)]
# save(sst_lm, file = "sst_lm.rda")
load("sst_lm.rda")
# sst_ARMA <-
# sst[!is.na(Mean_Ampl) & !is.na(asst) & sd != 0,
# {
# a <- arima(Mean_Ampl, xreg = asst,
# order = c(1, 0, 1), method = "ML")
# list(estimate = coef(a)[[4]],
# se = sqrt(diag(a$var.coef))[4])
# },
# by = .(lon, lat, month)]
# save(sst_ARMA, file = "sst_ARMA.rda")
sst_lm[, month := factor(month, levels = c(12, 1:11), ordered = T)]## lon lat month estimate se
## 1: 165.00 -76.070 1 3.465502 7.369890
## 2: 168.75 -76.070 1 6.279632 7.834251
## 3: 172.50 -76.070 1 9.603780 9.051243
## 4: 176.25 -76.070 1 13.800397 9.930859
## 5: 180.00 -76.070 1 7.110357 7.979328
## ---
## 35617: 333.75 87.159 12 -418.695999 589.880315
## 35618: 337.50 87.159 12 -66.482428 156.488025
## 35619: 341.25 87.159 12 -67.993298 205.546154
## 35620: 345.00 87.159 12 -288.471974 677.332404
## 35621: 356.25 87.159 12 -1405.224116 3258.036442sst_lm[, t := abs(estimate)/se]## lon lat month estimate se t
## 1: 165.00 -76.070 1 3.465502 7.369890 0.4702243
## 2: 168.75 -76.070 1 6.279632 7.834251 0.8015613
## 3: 172.50 -76.070 1 9.603780 9.051243 1.0610454
## 4: 176.25 -76.070 1 13.800397 9.930859 1.3896479
## 5: 180.00 -76.070 1 7.110357 7.979328 0.8910972
## ---
## 35617: 333.75 87.159 12 -418.695999 589.880315 0.7097982
## 35618: 337.50 87.159 12 -66.482428 156.488025 0.4248404
## 35619: 341.25 87.159 12 -67.993298 205.546154 0.3307933
## 35620: 345.00 87.159 12 -288.471974 677.332404 0.4258942
## 35621: 356.25 87.159 12 -1405.224116 3258.036442 0.4313101ggplot(sst_lm, aes(lon, lat)) +
geom_raster(aes(fill = estimate)) +
# geom_contour(aes(z = abs(estimate/se)), breaks = 2:100, color = "black") +
geom_point(size = 0.1, shape = 3, alpha = 0.4,
data = sst_lm[t > 2]) +
world3 + coord_quickmap() +
theme_elio +
facet_wrap(~month, labeller = labeller(month = month.abb), ncol = 3) +
scale_fill_gradient2(high = muted("red"), low = muted("blue"),
name = "Regresión (LM)", limits = c(-40, 40))
SST*
sst_comp <- sst[!is.na(Anom),
.(asst = mean(asst, na.rm = T)),
by = .(lat, lon, Anom, month)]
tmp <- sst_comp[, .(Anom = "Diferencia",
asst = asst[Anom == "Inactivo"] - asst[Anom == "Activo"]),
by = .(lat, lon, month)]
sst_comp <- rbind(sst_comp, tmp)
remove(tmp)
gdata <- sst_comp[Anom %in% c("Activo", "Inactivo")]
g <- ggplot(gdata, aes(lon, lat)) +
geom_contour(aes(z = asst, color = ..level..), binwidth = 2.5) +
world3 +
coord_quickmap() +
scale_color_gradient2(high = muted("red"), low = muted("blue")) +
theme_elioDEF: (líneas = 2.5)
g + facet_grid_paginate(month~Anom, labeller = labeller(month = month.abb),
nrow = 3, ncol = 2,
page = 1)
MAM
g + facet_grid_paginate(month~Anom, labeller = labeller(month = month.abb),
nrow = 3, ncol = 2,
page = 2)
JJA
g + facet_grid_paginate(month~Anom, labeller = labeller(month = month.abb),
nrow = 3, ncol = 2,
page = 3)
SON:
g + facet_grid_paginate(month~Anom, labeller = labeller(month = month.abb),
nrow = 3, ncol = 2,
page = 4)
No se puede ver a ojo mucha diferencia. Restando ambos campos (íneas = 0.5):
lineas <- seq(-3, 3, by = 0.5)
gdata <- sst_comp[Anom == "Diferencia"]
ggplot(gdata, aes(lon, lat)) +
geom_contour(aes(z = asst, color = ..level..), breaks = lineas) +
world3 +
coord_quickmap() +
facet_wrap(~month, labeller = labeller(month = month.abb), ncol = 3) +
scale_color_gradient2(high = muted("red"), low = muted("blue")) +
theme_elio
En algunos meses se ve que la zona del pacífico ecuatorial tiene grandes diferencias, con nomviembre, diciemre y enero caracterizado por diferencias positivas (años activos con mayor temperatura que inactivos) aunque lo contrario en septiembre y febrero.
Al hacer la clasificación binaria entre activos e inactivos se pierde información. Una alternativa es hacer la regresión de la SST en la amplitud de la onda QS3. Es decir, modelar \(QS3 = m\times SST + b + \epsilon\) para cada punto del globo.
invisible(
sst[!is.na(Mean_Ampl) & !is.na(asst), sd := sd(asst), by = .(lat, lon, month)]
)
# asst_lm <-
# sst[!is.na(Mean_Ampl) & !is.na(asst) & sd != 0,
# {
# a <- summary(lm(Mean_Ampl ~ asst))
# list(estimate = a$coefficients[2, 1],
# se = a$coefficients[2, 2])
# },
# by = .(lon, lat, month)]
# save(asst_lm, file = "asst_lm.rda")
load("asst_lm.rda")
# asst_ARMA <-
# sst[!is.na(Mean_Ampl) & !is.na(asst) & sd != 0,
# {
# a <- arima(Mean_Ampl, xreg = asst,
# order = c(1, 0, 1), method = "ML")
# list(estimate = coef(a)[[4]],
# se = sqrt(diag(a$var.coef))[4])
# },
# by = .(lon, lat, month)]
# save(asst_ARMA, file = "asst_ARMA.rda")
load("asst_ARMA.rda")
invisible({
sst_ARMA[, model := "ARMA"]
sst_lm[, model := "lm"]
})
sst_regr <- rbind(sst_ARMA, sst_lm)
remove(sst_ARMA, sst_lm)
invisible(
sst_regr[, t := abs(estimate/se)]
)
ggplot(sst_regr[model == "lm" & abs(estimate) < 40], aes(lon, lat)) +
geom_raster(aes(fill = estimate)) +
# geom_contour(aes(z = abs(estimate/se)), breaks = 2:100, color = "black") +
geom_point(size = 0.1, shape = 3, alpha = 0.4,
data = sst_regr[model == "lm" & t > 2]) +
world3 + coord_quickmap() +
theme_elio +
facet_wrap(~month, labeller = labeller(month = month.abb), ncol = 3) +
scale_fill_gradient2(high = muted("red"), low = muted("blue"),
name = "Regresión (LM)", limits = c(-40, 40))
Si quiero ser un poco más riguroso y tener en cuenta la estructura de correlación que seguro existe en las series, puedo usar un modelo ARMA(1, 1) para estimar la pendiente.
ggplot(sst_regr[model == "ARMA" & abs(estimate) < 40], aes(lon, lat)) +
geom_raster(aes(fill = estimate)) +
# geom_contour(aes(z = abs(estimate/se)), breaks = 2:100, color = "black") +
geom_point(size = 0.1, shape = 3, alpha = 0.4,
data = sst_regr[model == "ARMA" & t > 2]) +
world3 + coord_quickmap() +
theme_elio +
facet_wrap(~month, labeller = labeller(month = month.abb), ncol = 3) +
scale_fill_gradient2(high = muted("red"), low = muted("blue"),
name = "Regresión (ARMA)", limits = c(-40, 40))
En ambas figuras las cruces marcan regiones donde la magnitud de la regresión es mayor a 2 veces el error estandar de la misma (~ p-valor < 0.05).
(Alternativa: hacer test t-student para cada punto de latitud entre años activos e inactivos. )
# test <-
# sst[!is.na(Mean_Ampl) & !is.na(asst) & sd != 0,
# t.test(asst[Anom == "Activo"], y = asst[Anom == "Inactivo"]),
# by = .(lat, lon, month)]
#
# ggplot(test, aes(lon, lat)) +
# geom_raster(aes(fill = estimate)) +
# geom_point(aes(alpha = p.value < 0.05), size = 0.1, shape = 3) +
# scale_alpha_discrete(range = c(0, 0.7)) +
# scale_fill_gradient2(high = muted("red"), low = muted("blue")) +
# facet_wrap(~month, labeller = labeller(month = month.abb), ncol = 3) +
# coord_quickmap() +
# world3 +
# theme_elio